//
//  DEMConverter.cpp
//  dem-converter
//
//  Created by Giorgio Tino on 2/25/14.
//
//

#include "DEMConverter.h"

#include <libpng/png.h>
#include <stdlib.h>
#include <stdint.h>

#include <plib/ul.h>
#include <plib/ssgLocal.h>

#define FILE_ROW_STRING_LENGTH 80

#define isOutOfBounds(x)(x < 0) // Who knows if this will change in the future...

typedef struct
{
    int x;
    int y;
} iCoord;

typedef enum
{
    UNKNOWN = -1,
    NORTH,
    EAST,
    SOUTH,
    WEST
} eDirection;

DEMConverter::~DEMConverter()
{
    
}

bool DEMConverter::convert(const char *input, const char *output)
{
    FILE* f;
    char fileRow[FILE_ROW_STRING_LENGTH];
    char ch1,ch2;
    char id[FILE_ROW_STRING_LENGTH];
    int c = 0;
    int r = 0;
    int east = 0;
    int north = 0;
    float ah = 0;
    float step = 25.0;
    int units = 0; // METERS
    
    float* heights = NULL;
    
    /* LOADER */
    if((f = fopen(input,"r")) != NULL)
    {
		/* jump beyond header */
		fgets(fileRow,FILE_ROW_STRING_LENGTH,f);
		fgets(fileRow,FILE_ROW_STRING_LENGTH,f);
        sscanf(fileRow,"%s",id);
		fgets(fileRow,FILE_ROW_STRING_LENGTH,f);
		fgets(fileRow,FILE_ROW_STRING_LENGTH,f);
        
		/* read the number of columns */
		fgets(fileRow,FILE_ROW_STRING_LENGTH,f);
		sscanf(fileRow,"%c%c%d",&ch1,&ch2,&c);
        
		/* read the number of rows */
		fgets(fileRow,FILE_ROW_STRING_LENGTH,f);
		sscanf(fileRow,"%c%c%d",&ch1,&ch2,&r);
        
		/* read the x_start_coord (and  stores it in x)  */
		fgets(fileRow,FILE_ROW_STRING_LENGTH,f);
		sscanf(fileRow,"%d",&east);
        
		/* read the y_start_coord (and  stores it in y) */
		fgets(fileRow,FILE_ROW_STRING_LENGTH,f);
		sscanf(fileRow,"%d",&north);
        
		/* read the average height */
		fgets(fileRow,FILE_ROW_STRING_LENGTH,f);
		sscanf(fileRow,"%f",&ah);
        
		/* jump over another line */
		fgets(fileRow,FILE_ROW_STRING_LENGTH,f);
        
		/* read the lines containing the height values */
		/* This matrix is filled per columns, from the bottom to the top,
         and from the left to the right. So data in memory have the 'same'
         position (each other) than in the real world...*/
        
        /* allocate a height map of r * c */
        heights = new sgFloat[r*c];
        
        int i = 0;
		while(fgets(fileRow,FILE_ROW_STRING_LENGTH,f) != NULL)
        {
            sscanf(fileRow,"%f",&(heights[i]));
            i++;
        }
        
		/* close the file stream  */
		fclose(f);
        
        /* SAVE */
        FILE *fd = fopen ( output, "wb" ) ;

        if ( fd == NULL )
        {
            perror ( output ) ;
            ulSetError ( UL_FATAL,
                        "dem-converter: Failed to open '%s' for writing.", output ) ;
            return false ;
        }

        _ssgWriteString(fd, id);
        _ssgWriteInt( fd, r );
        _ssgWriteInt( fd, c );
        _ssgWriteFloat( fd, step );
        _ssgWriteInt( fd, units );
        _ssgWriteInt( fd, east);
        _ssgWriteInt( fd, north );

        //
        // Save height map
        //
        for(int i = 0; i < r*c; i++)
        {
            _ssgWriteFloat(fd, heights[i]);
        }
        
        delete[] heights;
        heights = 0;
        
        fclose ( fd ) ;
    }
	else
    {
		perror ( input ) ;
        ulSetError ( UL_FATAL,
                    "dem-converter: Failed to open '%s' for reading.\n", input ) ;
        return false ;
    }
    
    return true;
}

sgFloat** DEMConverter::load(const char* input, int* rows, int* columns, sgFloat* min, sgFloat* max)
{
    sgFloat** heightMap;
    
    char* id;
    int c = 0;
    int r = 0;
    int east = 0;
    int north = 0;
    int i = 0;
    int j = 0;
    int units = 0;
    float step;
    sgFloat hMin = MAXFLOAT;
    sgFloat hMax = -MAXFLOAT;
    
    /* LOADER */
    FILE *f = fopen ( input, "rb" ) ;
    
    if ( f == NULL )
    {
        perror ( input ) ;
        ulSetError ( UL_FATAL,
                    "load: Failed to open '%s' for reading.", input ) ;
        return NULL;
    }
    
    printf("\n\nloading : %s\n\n",input);
    
    _ssgReadString(f, &id);
    _ssgReadInt( f, &r );
    _ssgReadInt( f, &c );
    _ssgReadFloat( f, &step );
    _ssgReadInt( f, &units );
    _ssgReadInt( f, &east);
    _ssgReadInt( f, &north );
    
    /* allocate a height map of r * c */
    heightMap = new sgFloat*[r];
    heightMap[0] = new sgFloat[r*c];
    
    for(int i=1; i<r; i++)
    {
        heightMap[i] = &heightMap[0][i*c];
    }
    
    /* read the lines containing the height values */
    /* This matrix is filled per columns, from the bottom to the top,
     and from the left to the right. So data in memory have the 'same'
     position (each other) than in the real world...*/
    int n = r*c;
    while(n>0)
    {
        float value;
        if(j < r)
        {
            _ssgReadFloat(f,&value);
            heightMap[r-j-1][i] = value;
            j++;
        }
        else
        {
            j = 0;
            i += 1;
            _ssgReadFloat(f,&value);
            heightMap[r-j-1][i] = value;
            j++;
        }
        
        if (value > 0)
        {
            hMin = value < hMin ? value : hMin;
            hMax = value > hMax ? value : hMax;
        }
        
        n--;
    }
    
    /* close the file stream  */
    fclose(f);
    
    //print(heightMap, r, c);
    //printBox(heightMap, r, c);
    
    *rows = r;
    *columns = c;
    *min = hMin;
    *max = hMax;
    
    return heightMap;
}

int DEMConverter::saveAsPng(sgFloat** heightMap, int r, int c, sgFloat min, sgFloat max, const char* path) const
{
    FILE * fp;
    png_structp png_ptr = NULL;
    png_infop info_ptr = NULL;
    size_t x, y;
    png_byte ** row_pointers = NULL;
    /* "status" contains the return value of this function. At first
     it is set to a value which means 'failure'. When the routine
     has finished its work, it is set to a value which means
     'success'. */
    int status = -1;
    /* The following number is set by trial and error only. I cannot
     see where it it is documented in the libpng manual.
     */
    int pixel_size = 3;
    int depth = 8;
    
    fp = fopen (path, "wb");
    if (! fp) {
        goto fopen_failed;
    }
    
    png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
    if (png_ptr == NULL) {
        goto png_create_write_struct_failed;
    }
    
    info_ptr = png_create_info_struct (png_ptr);
    if (info_ptr == NULL) {
        goto png_create_info_struct_failed;
    }
    
    /* Set up error handling. */
    
    if (setjmp (png_jmpbuf (png_ptr))) {
        goto png_failure;
    }
    
    /* Set image attributes. */
    
    png_set_IHDR (png_ptr,
                  info_ptr,
                  c,
                  r,
                  depth,
                  PNG_COLOR_TYPE_RGB,
                  PNG_INTERLACE_NONE,
                  PNG_COMPRESSION_TYPE_DEFAULT,
                  PNG_FILTER_TYPE_DEFAULT);
    
    /* Initialize rows of PNG. */
    
    row_pointers = (png_byte **)png_malloc (png_ptr, r * sizeof (png_byte *));
    for (y = 0; y < r; ++y) {
        png_byte *row =
        (png_byte *)png_malloc (png_ptr, sizeof (uint8_t) * c * pixel_size);
        row_pointers[y] = row;
        for (x = 0; x < c; ++x) {
//            pixel_t * pixel = pixel_at (bitmap, x, y);
//            *row++ = pixel->red;
//            *row++ = pixel->green;
//            *row++ = pixel->blue;
            sgFloat height = heightMap[y][x];
            height = height > 0 ? height : 0;
            
            uint8_t value = (uint8_t)(( (height - min) / (max - min)) * 256);
            
            *row++ = value;
            *row++ = value;
            *row++ = value;
        }
    }
    
    /* Write the image data to "fp". */
    
    png_init_io (png_ptr, fp);
    png_set_rows (png_ptr, info_ptr, row_pointers);
    png_write_png (png_ptr, info_ptr, PNG_TRANSFORM_IDENTITY, NULL);
    
    /* The routine has successfully written the file, so we set
     "status" to a value which indicates success. */
    
    status = 0;
    
    for (y = 0; y < r; y++) {
        png_free (png_ptr, row_pointers[y]);
    }
    png_free (png_ptr, row_pointers);
    
    png_failure:
    png_create_info_struct_failed:
        png_destroy_write_struct (&png_ptr, &info_ptr);
    png_create_write_struct_failed:
        fclose (fp);
    fopen_failed:
        return status;
}

void DEMConverter::print(sgFloat** heightMap, int r, int c) const
{
    printf("\n");
    
    for(int i=0; i<r;i++)
    {
        for(int j=0; j<c;j++)
        {
            float height = heightMap[i][j];
            
            if(height < 0)
                printf("X ");
            else
                printf("%.0f ",height);
        }
        
        printf("\n");
    }
    
    printf("\n");
}

void DEMConverter::printVisited(bool** visited, int r, int c) const
{
    printf("\n");
    
    for(int i=0; i<r;i++)
    {
        for(int j=0; j<c;j++)
        {
            bool v = visited[i][j];
            
            if(v)
                printf("X ");
            else
                printf("O ");
        }
        
        printf("\n");
    }
    
    printf("\n");
}

void DEMConverter::printBox(sgFloat** heightMap, int r, int c) const
{
    int nVertices = (c+r)*4;
    
    if(nVertices < 6)
        return;
    
    printf("[%d %d] : %d\n", r, c, nVertices);
    
    iCoord current = {0 , 0};
    
    iCoord start = {-1, -1};
    
    //
    // Find the borders of the map to create the boarder box. A point is considered to be part of the border
    // if at least one of the adjacent points (top left bottom right) height is < 0 (DEM files have values of -9999 usually)
    //
    // NOTE : this is assuming that there is only one 'island' and no holes or crazy configurations! Let's keep it simple...
    //
    
    //
    // STEP 1 : find the strating point
    //
    bool found = false;
    
    for(int i = 0; i < r; i++)
    {
        for(int j = 0; j < c; j++)
        {
            printf("i,j is [%d,%d] is %f : %d\n",i,j,heightMap[i][j],isOutOfBounds(heightMap[i][j]));
        }
    }
    
    for(int i = 0; (i < r) && !found; i++)
    {
        for(int j = 0; (j < c) && !found; j++)
        {
            printf("i,j is [%d,%d] is %f : %d\n",i,j,heightMap[i][j],isOutOfBounds(heightMap[i][j]));
            
            if(!isOutOfBounds(heightMap[i][j]))
            {
                start.x = j;
                start.y = i;
                found = true;
            }
        }
    }
    
    //
    // Completely empty map. Very weird, should NEVER happen!
    //
    if(!found)
    {
        printf("NOT FOUND!\n");
        return;
    }
    
    bool** visited = new bool*[r];
    visited[0] = new bool[r*c]();
    
    for(int i=1; i<r; i++)
    {
        visited[i] = &visited[0][i*c];
    }
    
    printf("Start : [%d,%d]\n",start.x,start.y);
    
    visited[start.y][start.x] = true;
    
    //
    // STEP 1 : from the starting point, start following the boundary (top left bottom right) always excluding the direction we are coming from to avoid loops. When the current position is equal to the starting position, we consider the box closed. NOTE : again, this works for convex trrains with one island only and no holes. Which is correct for all the DEM files I have ever seen.
    //
    current = start;
    
    eDirection direction = UNKNOWN;
    
    do
    {
        iCoord next = current;
        
        //printf("{%d,%d}\n",next.x,next.y);
        
        // Check NORTH
        if((direction != SOUTH) &&
           (current.y-1 >= 0) &&
           _isBorder(heightMap, visited, r, c, current.x , current.y-1))
        {
            next.x = current.x;
            next.y = current.y - 1;
            
            // next is north!
            direction = NORTH;
        }
        else if((direction != WEST) &&
                (current.x+1 < c) &&
                _isBorder(heightMap, visited, r, c, current.x+1, current.y))
        {
            // next is east!
            next.x = current.x + 1;
            next.y = current.y;
            
            direction = EAST;
        }
        else if((direction != NORTH) &&
                (current.y+1 < r) &&
                _isBorder(heightMap, visited, r, c, current.x, current.y+1))
        {
            // next is south!
            next.x = current.x;
            next.y = current.y+1;
            
            direction = SOUTH;
        }
        else if((direction != EAST) &&
                (current.x-1 >= 0) &&
                _isBorder(heightMap, visited, r, c, current.x-1, current.y))
        {
            // next is west!
            next.x = current.x - 1;
            next.y = current.y;
            
            direction = WEST;
        }
        else
        {
            direction = UNKNOWN;
        }
        
        if(direction != UNKNOWN) // should never happen
        {
            printf("[%d,%d]\n",next.x,next.y);
            visited[next.y][next.x] = true;
            
            printVisited(visited, r, c);
        }
        else
        {
            break;
        }
        
        current = next;
    }
    while(current.x != start.x || current.y != start.y);
}

bool DEMConverter::_isBorder(sgFloat** heightMap, bool** visited, int r, int c, int x, int y) const
{
    if(isOutOfBounds(heightMap[y][x]) || visited[y][x])
        return false;
    
    // Top left
    for(int i = y-1; i < y+2; i++ )
    {
        for(int j = x-1; j < x+2; j++)
        {
            if((i >= 0) &&
               (j >= 0) &&
               (i < r) &&
               (j < c) &&
               (i != y || j != x) &&
               isOutOfBounds(heightMap[i][j]) &&
               !visited[i][j])
            {
                return true;
            }
        }
    }
    
    return false;
}